Three-dimensional vortex structure of a fast rotating Bose-Einstein condensate with 

harmonic-plus-quartic confinement 
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We address the challenging proposition of using real experimental parameters in a three- 
dimensional (3D) numerical simulation of fast rotating Bose-Einstein condensates. We simulate 
recent experiments [V. Bretin, S. Stock, Y. Seurin and J. Dalibard, Phys. Rev. Lett. 92, 050403 
(2004); S. Stock, V. Bretin, S. Stock, F. Chevy and J. Dalibard, Europhys. Lett. 65, 594 (2004)] 
using an anharmonic (quadratic-plus-quartic) confining potential to reach rotation frequencies (Q) 
above the trap frequency (ui). Our numerical results are obtained by propagating the 3D Gross- 
Pitaevskii equation in imaginary time. For < ui± , we obtain an equilibrium vortex lattice similar 
(as the size and number of vortices) to experimental observations. For > lu±_ we observe the 
evolution of the vortex lattice into an array of vortices with a central hole. Since this evolution 
was not visible in experiments, we investigate the 3D structure of vortex configurations and 3D 
effects on vortex contrast. Numerical data are also compared to recent theory [D. E. Sheehy and 
L. Radzihovsky, Phys. Rev. A 70, 063620 (2004)] describing vortex lattice inhomogeneities and a 
remarkably good agreement is found. 

PACS numbers: 03.75.Fi,02.70.-c 



I. INTRODUCTION 

In recent years, several experimental studies provided 
evidence for the existence of quantized vortices in rotat- 
ing Bose-Einstein condensates (becs) UBHSSIE0- 

The condensate is typically confined by a harmonic 
(quadratic) potential with transverse frequency u>±_ and 
starts to nucleate vortices when the rotation frequency O 
exceeds a critical value Vt c . For increasing Q > f2 c , more 
and more vortices appear and arrange themselves into a 
regular triangular (Abrikosov) lattice. 

The fast rotation regime, corresponding to O > u>±, is 
particularly interesting to explore since a rich variety of 
scenarios are theoretically predicted: formation of giant 
(multi-quantum) vortices, vortex lattice melting or quan- 
tum Hall effects. This regime is experimentally delicate 
to investigate [j| since for = uo±_ the centrifugal force 
compensates the trapping force and the confinement van- 
ishes. Using evaporative spin up, the Boulder group has 
recently created condensates with rotation frequencies of 
the order of 0.99wj_ and studied the properties of the 
vortex lattice in the lowest Landau level 0, El El • 

Another experimental approach to reach the fast 
rotation regime was explored by the Ecole Normale 
Superieure (ENS) group 0, [lJ It consists in 

modifying the quadratic trapping potential by super- 
imposing a blue detuned laser beam to the magnetic 
trap holding the atoms. The resulting harmonic-plus- 
Gaussian potential removes the singularity at the limit 
fl = lu±_ and allows to reach rotation rates up to Q ~ 
1.05loj_. The trapping potential used in these exper- 
iments can be well approximated by a quadratic-plus- 
quartic form, which has been extensively studied lately 
Ellllllllllllllllllllilllllii]. Different transi- 
tions involving a rich variety of vortex states are theoret- 
ically predicted when f2 is increased: from a dense vortex 



lattice to an array of singly quantized vortices with a cen- 
tral hole and, finally, to a giant (multiquantum) vortex 
or directly from a vortex lattice to a giant vortex. 

For the highest rotation rates reached in experiments, 
neither giant vortices nor vortex arrays with hole were 
clearly observed 0,0]. In exchange, a dramatic change 
in the appearance of the condensate was reported: the 
vortices are less visible even thought the gas remains ul- 
tracold and in fast rotation. The most likely explanation 
for this intriguing behavior was related to the transient 
character of the observed states leading to a fragile vortex 
lattice dominated by three-dimensional (3D) effects (vor- 
tices appear to have some excitation or bending leading 
to poor optical contrast). 

Since such effects are not trackable with previous (2D) 
numerical approaches, the purpose of the present contri- 
bution is to investigate the 3D structure of such conden- 
sates by numerically generate the corresponding Gross- 
Pitaevskii (GP) wave function. This is not without 
its challenges, since the description of a prolate (cigar- 
shaped) condensate with a large number of vortices (ex- 
ceeding 100) requires very high spatial resolution and ac- 
curate integration schemes. Computations become very 
expensive at high rotation frequencies, which explains 
why such 3D simulations are not, to the author's knowl- 
edge, currently available in the open literature. 

The numerically generated 3D condensates can be seen 
in Fig. n For increasing rotation frequencies, the vortex 
lattice evolves to a vortex array with a hole, which con- 
firms the scenario theoretically predicted JTa H H HJ 
and also observed in 2D simulations [TsL l24j |. Since such 
transition was not observed in experiments, we qualita- 
tively analyze the obtained vortex states, with a par- 
ticular emphasize on the 3D features of vortex merging 
leading to a central hole in the condensate. 

Our analysis is then extended to quantitative compar- 
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isons to experiments and theoretical predictions. We first 
check that physical parameters (size, chemical potential) 
of numerical condensates correspond well to available ex- 
perimental ones. We show in particular that the rotation 
frequencies reached in experiments were not enough high 
to obtain an annular condensate. We also measure from 
simulations the intervortex spacing and compare the nu- 
merical results to recent theory of Sheehy and Radzi- 
hovsky [2(|II3 describing vortex lattice inhomogeneities. 
A remarkably good agreement is found. Finally, we dis- 
cuss how 3D structure of vortices can affect optical con- 
trast of transient states observed in experiments. 



II. PHYSICAL PARAMETERS AND 
NUMERICAL APPROACH 

We consider a BEC of N atoms confined by the trap- 
ping potential V and rotating along the z axis at angu- 
lar velocity ft. In the experiments at ENS [T^. ll3L Il4| . 
N = 3xl0 5 atoms and the trapping potential can be 
written as the superposition of the harmonic potential 
Vh created by the magnetic trap and the potential U(r) 
introduced by the laser beam propagating along the z 



V(r,z) = V h (r,z) + U(r), 



(1) 



with r = \J x 2 + y 2 and 



1 



1 



K = -m(cf )V + im^V, C/(r) = C/ e- 2 ^ . 

(2) 

The trapping frequencies are a;? — 2nx75.5 Hz and 
lu z = 27rxll Hz, resulting in a cigar-shaped condensate. 
The laser waist is w — 25 /xm and the amplitude of the 
laser beam is Uq = fcsx90 nK. 

For r/w sufficiently small, the potential V(r) can be 
approximated by: 



Vi 



1 1 (0)x2 2U 



2U Q 



w 1 



\ muj2 z z 2 . (3) 



For this quadratic-plus- quartic potential, the transverse 
trapping frequency is decreased to u>± — 27rx65.6 Hz. 
Since the amplitude Uq of the laser beam is low in exper- 
iments, the quadratic part of the potential V\ remains 
positive (repulsive interactions) and the quartic part is 
very small. It is interesting to note that a stronger ampli- 
tude U could generate a quartic-minus- quadratic poten- 
tial, which was theoretically studied in Refs. [2fll22 . l2fij . 

The numerical results presented in this paper were 
obtained using a quadratic-plus- quartic potential (Eq. 
EDl, for which extensive theoretical studies are available 
EE EE E3, EE IH, E3, E3, ISE E3| Numerical sim- 
ulations using the quadratic-plus-Gaussian original po- 
tential (Eq. |5J) showed the same qualitative evolution of 
the vortex configuration as in Fig. ^ with a transition 
to a vortex array with hole for a slightly lower rotation 
frequency. 



As a numerical approach, we compute the macro- 
scopic wave-function ip(x, y, z) by propagating the three- 
dimensional Gross-Pitaevskii (GP) equation in i maginary 
time by the numerical method used in Refs. 122. 123 . l29 | . 
After rescaling the GP equation as in Ref. 30], a hy- 
brid Runge-Kutta-Crank-Nicolson scheme is used for the 
time integration and a sixth-order compact finite differ- 
ence scheme for the space discretization. 

As initial condition we generally use a vortex-free den- 
sity distribution following the Thomas-Fermi (TF) law: 



p TF (r, z) 



4irh 2 c 



fJ ,-V 1 (r,z) + -mQ 2 r 2 ) , (4) 



where a s — 5.2 nm is the scattering length and \x the 
chemical potential given by the constraint J d 3 rp TF = N. 
For the quadratic-plus- quartic trapping potential V±, an 
exact analytical form for /i can be derived I5j depending 
on the value of f2 which dictates the shape of the conden- 
sate (with or without a hole). The maximum transverse 
radius R± and longitudinal half-length R z can be then 
calculated from Eq. m order to estimate the dimen- 
sions of the rectangular computational domain. For high 
f2 (when the condenstae is nearly spherical and more than 
100 vortices are present), up to 240 x 240 x 240 grid points 
are used to compute equilibrium states. 

The post-proces sing of the results follows the exper- 
imental approach |lll (with the difference that the 
radial expansion after the time of flight is not numer- 
ically simulated). The numerical 3D wave- function is 
converted to an atomic density p(x,y,z) = \ip(x, y, z)\ 2 
and integrated along the rotation (z) axis. The result- 
ing 2D-density p z (x, y) (isocontours are displayed in Fig. 
□ last row of images) will be used in the following for 
comparison to experiments and theory. 



III. DESCRIPTION OF THE RESULTS 

The evolution of the 3D structure of the condensate 
with increasing Q, can be seen in Fig. ^ We start with 
a qualitative description of vortex configurations. The 
obtained results will be then analyzed quantitatively and 
compared to available experimental and theoretical val- 
ues. All quantitative parameters discussed in this paper 
are summarized in table [I] 



A. Vortex configurations 

For rotation frequencies below ui± (fi/(27r) = 60 and 
64) the condensate has the usual prolate shape (see Fig. 
^ first two columns). Vortices near the center of the con- 
densate are straight and form a regular triangular lattice. 
Vortices located near r — R± are bending, reaching the 
surface of the condensate using the shortest path. These 
outer vortices are not symmetrically arranged and have 
different lengths. It is interesting to note that for these 
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FIG. 1: (Color online) Numerically generated condensates obtained using a quadratic+quartic trapping potential with the 
parameters corresponding to experiments of fl2T.ll3jl. Each column corresponds to a value of the rotation frequency - from left 
to right: 0,/2-k = 60,64,66,70.6,73 (respectively, Q/uj± = 0.92,0.98,1.01,1.08,1.11). Three-dimensional views of the vortex 
lattice identified by means of iso-surfaces of low atomic-density (first two rows) and contours of density integrated along the 
rotation (z) axis. Note that the formation of the hole in the condensate is not complete for Q/2tt — 70.6 and we still distinguish 
individual singly-quantized vortices in the center (see also Fig. [5] for a detailed picture of this configuration). 



Q/(2n) 


60 
0.92 


64 
0.98 


66 
1.01 


70.6 
1.08 


73 
1.11 


R± \jj/m] 


10.4 


12.2 


13.2 


17.2 


19.2 


R z [fim] 


29.0 


25.4 


22.5 


20.1 


18.6 


N v 


37 


51 


62 


126 


113 


L z [units of h] 


17.4 


28.5 


39.1 


122.6 


239.1 




2.15 


1.84 


1.65 


1.36 


1.76 



TABLE I: Summary of the characteristics of numerically gen- 
erated condensates: (maximum) transverse radius R± and 
longitudinal half-length R z ; number of vortices N v and an- 
gular momentum L z = i J d s rip (ydip/dx — xdip/dy); scaling 
constant for the ratio between vortex-core radius r v and heal- 
ing length £ [obtained from integrated density p z (x,y)]. 



two values of f2, the number of vortices N v we find nu- 
merically (N v = 37 and 51) is very close to that visible 
in experimental pictures [lj (N% xpt = 30 and 52). 
Starting with Sl/(2n) = 66 (fl/oj± = 1.01), the exper- 



imental pictures show a lack of contrast for entire zones 
of the vortex lattice. Vortices are less visible and do not 
allow a proper estimation of the rotation frequency from 
vortex surface density. Numerical condensate for this 
rotation frequency (Fig. ^ third column) display a well- 
defined triangular vortex lattice. Most of the vortices are 
straight and join the top and bottom ends of the conden- 
sate which are almost flat. This particular shape of the 
condensate corresponds well to that predicted from the 
TF law (@J. Indeed, for ft = uj±, the density distribution 
p TF (r, z) depends only on the quartic part of the trapping 
potential V\ and the surface of the condensate defined as 
{Ptf = 0} is flat near the rotation z axis. 

For rotation frequencies exceeding cjj^, experimental 
condensate exhibits a local minimum in the central den- 
sity, but the theoretically predicted 0, 13 transition to 
a vortex lattice with a hole (annular condensate) is not 
experimentally reported. This is the case in our simu- 
lations (Fig. The rotation frequency corresponding 
to this transition is found to be 0/j/(27r) = 71, a value 
close to the TF prediction ft£ F /(27r) = 70. These val- 
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ues are already larger than those attained in experiments 
(f2/(27r) < 69), which can simply explain why the hole 
was not experimentally observed. 

The numerically generated condensates before and af- 
ter transition to an annular condensate are shown in Fig. 
E] (last two columns of images). For fi/(27r) = 70.6, the 
central hole is not yet formed since the top and bottom 
depletions have not merged. At the center of the con- 
densate, the density is very low but not zero, and we 
can still distinguish individual vortices from isocontours 
of the density integrated along the z axis (Fig. Since 
the contrast in this last image is low near the center, 
we show details of the vortices near the rotation axis in 
Fig. [21 In the center there are three vortices with larger 
cores that start to reconnect at the top and bottom of 
the condensate. This merging process is highly three- 
dimensional and will finally lead to the formation of a 
central hole for higher Vt. 




FIG. 2: (Color online) Details of the vortex configuration 
for fi/(27r) = 70.6. Vortices near the rotation axis are iso- 
lated, showing the merging process that will finally lead to 
the formation of a central hole. Insert: top view of the same 
configuration. 

The structure of the condensate is completely different 
for Vl/(2tt) — 73 (last column of images in Fig. 2J. The 
condensate is nearly spherical, with a large central hole 
surrounded by three concentric circles of singly quan- 
tized vortices. Most of the 113 identified vortices are 
bent, reaching either inner or outer faces of the conden- 
sate. Since convergence for this case is particularly slow 
(two weeks of computational time is necessary using a 
PC workstation), we did not explore cases for higher f2. 
For the considered parameters, a second transition to a 
configuration with a pure giant vortex (without singly 
quantized vortices in the annular region) may occur at 
very high rotation frequencies |24j that are not numeri- 
cally affordable in 3D. 



B. Vortex lattice inhomogeneity 

We now turn on more quantitative analysis of numer- 
ical results. Before analyzing the characteristics of the 
vortex lattice, we first check that the dimensions of the 
numerical generated condensates correspond well to ex- 
perimental ones. The density p z is integrated along the 
azimuthal direction 9 to get the radial density profile 
p z ' e (r). This profile is fitted to the Thomas- Fermi dis- 
tribution (£Q), taking the chemical potential \i and the 
rotation frequency Q as adjustable parameters. The the- 
ory fit value of f2 is within 1% of the value of f2 for which 
the computation was done. 




QI2n 

FIG. 3: (Color online) Chemical potential (p) and maximum 
transverse radius of the condensate (R±) as functions of the 
rotation frequency. Experimental measurements from Ref. 
[T3 | (squares), numerical results (circles) and Thomas- Fermi 
theoretical prediction (solid line). 

The resulting chemical potential fj, and the transverse 
radius R± (which is the maximum radius for the con- 
densate with hole) are compared in Fig. |3|to experimen- 
tal values from Ref. [14j and Thomas-Fermi approxima- 
tion For the experimentally available range of rota- 
tion frequencies, numerical results are in good agreement 
with experimental and theoretical values. For values of 
fl not available experimentally, numerical results follow 
the TF prediction. In particular, the numerical value 
Uh/(2Tr) = 71 for which the central hole first appears 
in the condensate (corresponding to a chemical potential 
fj, = 0) is well predicted by the TF law (S1^ f /(2tt) = 70). 

We continue our dimensional analysis by extracting the 
characteristics of the vortex lattice: namely the intervor- 
tex spacing b v and the vortex core size r v . We follow a 
similar post-processing procedure as in Ref. Using 
the integrated (along z) density field p z (r, 9), we identify 
vortex centers by 2D searching of local minima. Result- 
ing points are checked to correspond to vortices visible in 
Fig. n (last row of images) . Assuming a triangular lattice 
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structure, we select vortices for which the six nearest vor- 
tex neighbors form a hexagonal pattern. Only for such 
vortices (i.e. vortices close to R± are discarded), is the in- 
tervortex spacing b v measured by averaging the distance 
from the vortex center to the centers of the six neighbors. 
The vortex core radius r v is measured as follows: for a 
given vortex located at (tq, 9q), the density profile p*(r) 
along the radius passing through the center of the vortex 
is extracted from the 2D field p z ; by subtracting the inte- 
grated TF density profile p^ F (r) [corresponding to Eq. 0] 
integrated along z] , we obtain a vortex-core residual that 
is fitted with a Gaussian profile: 



(n/27i=60) 



PtfM - P Z v( r ) = Aex P 



1 



(r - r Q ) 2 /r v < 



(5) 



The amplitude A is used to define the vortex contrast 
[Tl| as A/p* w (ro), i.e. the ratio between the "missing" 
column density at vortex center r and the corresponding 
TF value. Only vortices with a contrast greater than 0.7 
are considered to compute core radii r v . 

Figure 0] shows the variation of r v and b v as functions 
of the non-dimensional radius r/Rj_. Values are given in 
/im and rotation frequencies f2/wj_ < 1.01 are considered 
(condensates without hole). As expected 0, 0, the 
core radius r v scales with healing length, defined from the 
TF density fit £(r) = [Sna.p^r)}- 1 / 2 . The scaling con- 
stant (also summarized in Tab. P) decreases with f2, with 
values comparable to those found in Ref. 11] for a har- 
monic trapping potential. We recall that the values pre- 
sented here correspond to a post-processing for r v using 
integrated density p 2 , as in experiments. A similar post- 
processing using the 2D density field p extracted from 
the 3D simulation at z = 0, revealed scaling constants for 
r„/£ of order of 1 (more precisely, r„/£ ~ 0.98, 0.93, 0.86 
for, respectively, Q/2tt = 60,64,66). 

The calculated intervortex spacing b v is compared in 
Fig . H] to recent theory of Sheehy and Radzihovsky 
[26l l27j . They expressed the vortex density n v (r) as a 
function of the local superfluid density p s (r) : 
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FIG. 4: (Color online) Variation of vortex core radius r v and 
intervortex spacing b v (values in ^m) as functions of the non- 
dimensional radius r/R±. For each plot, the value of the 
rotation frequency (Q/2tt) is indicated. In plots displaying 
r v , solid line represents the variation of the healing length £, 
scaled by a constant indicated in the legend. Variation of b v 
is c omp ared to theory prediction of Sheehy and Radzihovsky 
HH H9 (solid line) and the estimation assuming a uniform 
(rigid-body) vortex distribution (dashed line). 



n v (r) 



+ In[h/(2.718 mQ£,t)] V 2 [ln{p s (r))] . (6) 



The second term in (JHJl is a small correction to the vortex 
density for a uniform vortex distribution corresponding 
to a rigid-body rotation n V Q = (flm)/(nh). The vortex 
density n v can be converted to intervortex spacing by: 



Mr) = ^2/(3^n v (r)). 



(7) 



Numerical results are compared to theoretical predic- 
tions using in Eq. [fj] the TF fit for the integrated den- 
sity profile (p s (r) = p^ F (r)) and the characteristic length 
for the vortex-core defined as [2(| = h/(mu±R±), 
The agreement is remarkably good. For S1/27T < 64, the 
density profile is close to an inverted parabola (the influ- 
ence of the quartic term being small) and b v is monoton- 
ically increasing with r. Similar results were reported 



for a harmonic trapping potential |ll| . As expected, 
the estimation using the rigid-body rotation assumption 
(dashed line in the plot) becomes better with increas- 
ing 11 (the lattice becomes denser). For H,/2tt — 70.6, 
the density profile has a Mexican-hat structure and vor- 
tices are constrained to agglomerate towards the center, 
where density is small. The intervortex spacing is small 
near the center and increases to the rigid-body value near 
r/R± ~ 0.5 where the density is maximum. The theory 
nicely illustrates this complex dependance of b v on the 
radial position. 



IV. DISCUSSION AND CONCLUSION 

We have presented in this paper three-dimensional 
numerical results for a fast-rotating BEC trapped in 
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quadratic-plus-quartic potential corresponding to exper- 
iments at ENS 0, The obtained vortex configu- 
rations show a transition from a dense vortex lattice to 
a vortex array with a central hole at a critical rotation 
frequency f2/j/(27r) = 71. This result confirms theoret- 
ical and 2D numerical results 0, EB E3 and goes be- 
yond experimental observations, since experiments failed 
to reach rotation frequencies close to Qh- 

Our results also support the assumption that vor- 
tices are less visible in experiments for Q/(2tt) > 66 be- 
cause of the fragility of the vortex lattice which becomes 
dominated by 3D-effects, such as vortex bending. In or- 
der to illustrate this statement it is worth describing how 
the condensate evolves in the "imaginary" time (i.e. how 
it relaxes to an equilibrium state). 

The imaginary-time evolution of the condensate looks 
similar to a real-time evolution. When suddenly increas- 
ing fl, new vortices are generated at the border of the con- 
densate and enter the condensate. In the first stages of 
the computation, 3D vortex lines are strongly distorted, 
giving a spaghetti image of the lattice. Close to equilib- 
rium, vortices become straight in their central part and 
arrange themselves in a more and more regular lattice. 
Convergence is particularly slow at the end of the com- 
putation when the position and shape of vortices evolve 
very slowly. Convergence is considered when the energy 
remains constant (relative fluctuations less than 10 -6 ) for 
a relatively long time to be sure that a stable state was 
obtained. The convergence time is much longer (roughly 
by a factor of 2) for values of the rotation frequency ex- 
ceeding u±. 




An example of intermediate states of the condensate 
before reaching a converged equilibrium state is displayed 
in Fig. |SJ The simulation corresponds to a quadratic- 
plus-Gaussian trapping potential PJl. (closer to the ex- 
perimental one) and a vortex configuration without hole. 
Transient states look closer to experimental pictures than 
the equilibrium states presented in Fig. ^ Three- 
dimensional exploration of the condensate reveals that 
vortices which are less visible have distorted structures 
which diminish the contrast in an integrated view along 
the z axis. These effects are stronger for condensates dis- 
playing a central depletion; even for equilibrium states 
of such condensates, it is difficult to distinguish individ- 
ual vortices in the center, as can be seen in Fig. 2] for 
f2/(27r) = 70.6. This confirms the hypothesis [l^ of the 
fragility of the experimental vortex lattice for high ro- 
tation frequencies: for transient states, 3D vortex lines 
have some excitations leading to a poor optical contrast. 
It is possible that the very low temperature in experi- 
ments slows down the dissipative process allowing only 
the observation of transient states dominated by 3D ef- 
fects. But is not to be excluded that a thermal excitation 
may be at the origin of the vortex-line bending responsi- 
ble for low optical contrast and, therefore, increasing the 
temperature in experiments is not a solution to improve 
vortex lattice contrast. 



Our simulations also offer a detailed 3D picture of vor- 
tex configurations that is not available from experiments 
and 2D simulations. In particular, the vortex merging 
leading to the formation of the central hole in a conden- 
sate is proved to be highly three-dimensional. Quantita- 
tive measurements of the intervortex spacing give a new 
validation of the theoretical study of Sheehy and Radz- 
ihovsky [2(| |53 predicting vortex lattice inhomogeneity 
from local density profile. An interesting question re- 
maining for future numerical investigations is whether or 
not the condensate trapped in a quadratic-plus-quartic 
potential enters the lowest Landau level regime for ~ 



FIG. 5: Example of energy decrease during the propagation of 
3D Gross-Pitaevskii equation in imaginary time. Simulation 
for f2/(27r) = 66, using the quadratic-plus-Gaussian trapping 
potential Energy is normalized by the equilibrium (final) 
value Ef. Inserts show iso-contours of the integrated (along 
z) density corresponding to three successive time instants rep- 
resented on the energy curve. 
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